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Abstract 


The radio-emitting neutron star population encompasses objects with 
spin periods ranging from milliseconds to tens of seconds. As they 
age and spin more slowly, their radio emission is expected to cease. 
We present the discovery of an ultra-long period radio-emitting neu- 
tron star, PSRJ0901—4046, with spin properties distinct from the 
known spin and magnetic-decay powered neutron stars. With a spin- 
period of 75.88 s, a characteristic age of 5.38 Myr, and a narrow 
pulse duty-cycle, it is uncertain how radio emission is generated and 
challenges our current understanding of how these systems evolve. 
The radio emission has unique spectro-temporal properties such as 
quasi-periodicity and partial nulling that provide important clues to 
the emission mechanism. Detecting similar sources is observationally 
challenging, which implies a larger undetected population. Our dis- 
covery establishes the existence of ultra-long period neutron stars, 
suggesting a possible connection to the evolution of highly magnetized 
neutron stars, ultra-long period magnetars, and fast radio bursts. 


1 Introduction 


Radio pulsars are rotation-powered neutron stars, which emit coherent beams 
of radio emission generated by highly relativistic particles in regions above 
their magnetic poles. Their known spin periods (P) range from 1.4ms to 
23.5s and they are divided into various sub-classes (e.g. rotating radio tran- 
sients, millisecond pulsars, magnetars - https://www.atnf.csiro.au/research/ 
pulsar/psrcat/) depending on their observational properties. Particle acceler- 
ation and abundant electron—positron pair production is postulated to be an 
essential condition for the coherent radio emission from pulsars, with the par- 
ticle acceleration potential expected to be lower for larger spin periods. As 
seen in most neutron stars, the radio emission is also expected to be strongly 
inhibited, or cease if the magnetic field configuration and strength, exceed 
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the quantum critical field (By. = 4.413 x 10'3G) [1]. Here, we present the 
discovery of a highly-magnetized 75.88s period radio-emitting neutron star, 
PSR J0901—4046 which challenges these conditions for, and the nature of, the 
radio emission and raises questions about the spin evolution of neutron stars 
in general. 


2 The Data 


2.1 The discovery and properties of PSR J0901—4046 


PSR J0901—4046 was a serendipitous single pulse discovery at 1284 MHz 
on 27 September 2020, in an observation directed at the high mass X-ray 
binary, Vela X-1, during simultaneous image and time domain searches by the 
Meer(more) TRAnsients and Pulsars (MeerTRAP - https://www.meertrap. 
org/) and ThunderKAT (http://www.thunderkat.uct.ac.za) projects at the 
MeerKAT radio telescope in South Africa. The pulse was initially detected 
in the MeerTRAP beamformed data in a single coherent tied-array beam of 
angular diameter ~ 45 arcseconds. A review of the MeerTRAP data for that 
observation revealed that there were further wide, but weaker pulses, which 
were missed by the real-time single pulse detection system. A total of fourteen 
pulses were identified in the beamformed time domain searches, which were 
regularly spaced in a span of ~ 30 minutes. A periodicity analysis resulted in 
an initial period of P = 75.89+0.01 seconds. The corresponding full time and 
frequency integration image of the field revealed an associated point source at 
the location of the coherent beam. These data were re-imaged at the smallest 
possible integration time of 8 seconds and more pulses were identified. An ini- 
tial inspection of the 8-second images from two other epochs where MeerTRAP 
data were not available, also revealed that the source exhibited a consistent 
periodicity. These snapshot images allowed the source to be localized to arcsec- 
ond precision. The deepest image of the field shows a partially visible, diffuse 
shell-like structure surrounding PSR J0901—4046 , which is possibly the super- 
nova remnant from the event that formed the neutron star. The complexity 
of the field in terms of diffuse emission requires additional analysis to deter- 
mine a robust association of this radio shell with PSR J0901—4046. No known 
pulsars are located within 2 degrees of this sky location. 

A total of six L-band (856 — 1712 MHz) and one UHF-band (544 — 1088 
MHz) observations have been performed between September 2020 and May 
2021. During these, we detect single pulses from every rotation of the source. 
The L-band data have resulted in the timing solution shown in Table 1. Despite 
the large jitter in the pulse shapes of single pulses, we obtain remarkably sta- 
ble pulse profiles over the various epochs due to the high signal-to-noise ratios 
(S/Ns). Using 29 times of arrival (ToAs), typically two per epoch, over 7.4 
months, we measure timing residuals with a low root-mean-square (rms) of 
5.7 ms (see Extended Data Figure 1). When compared to the pulse period, 
the fractional accuracy of ~ 7 x 10~° is comparable to the most accurately 
timed millisecond pulsars. We do not find any evidence of timing noise or 
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covariance of spin parameters with position. PSR JO901—4046 has a best fit 
dispersion measure (DM) of 52 + 1 pecm~® and average half-power pulse 
widths of ~ 300ms at both L- and UHF-band suggesting no evidence for 
radius-to-frequency mapping. We measure pulse-averaged peak flux densities 
of 89.3 + 2.7 mJy beam~+ and 169.3 + 14 mJy beam~! at L-band and UHF- 
band, respectively, with a period-averaged flux density of 408 +5 wJy beam + 
at L-band. The measured DM corresponds to distances of approximately 0.3 
and 0.5 kpc according to the YMw16 [2] and NE2001 [3] Galactic electron 
density models, respectively. The period (P = 75.88 s) and period derivative 
(P = 2.25 x 107!8 s s~1; pulsar spin-down rate) correspond to a characteris- 
tic age, surface magnetic field strength, and spin-down luminosity of 5.3 Myr, 
1.3 x 10!4 G and 2.0 x 1078 erg s~! assuming a dipolar magnetic field config- 
uration, respectively (see Figure 1). This discovery confirms the existence of 
ultra-long period neutron stars. 


3 Results 


3.1 Radio emission properties 


Single pulse analyses of the radio emission from PSRJ0901—4046 reveal 
remarkable and unusual spectro-tempo-polarimetric properties, quite unlike 
anything seen in known radio pulsars. We notice that the pulse shape is variable 
both inter-epoch and intra-epoch, but some features persist. Overall, the sin- 
gle pulses studied over 6 epochs can be grouped into 7 different types, namely: 
normal, quasi-periodic, spiky, double-peaked, partially nulling, split-peak and 
triple-peaked as shown in Figure 2. Although magnetars are sometimes seen 
to emit wide, bright radio pulses that comprise several sub-pulse components 
of varying widths and amplitudes, these are more chaotic within and between 
subsequent pulses. 

In some of the bright pulses we measure a quasi-periodicity in the sub-pulse 
components which at times appear to be harmonically related between pulses 
(see Extended Data Figure 2). In some others we see multiple quasi-periods 
within a single rotation as seen in Extended Data Figure 3. Overall, the quasi- 
periods are common across the UHF- and L-band observations. We observe 
the width of the sub-pulse components in PSR JO901—4046 to be exactly half 
of the quasi-period. The shortest and longest quasi-periods we measure are 
9.57 ms (104 Hz) and 338 ms (2.96 Hz) respectively (see Extended Data Figure 
4). Similar quasi-periodic features have been observed in fast radio bursts 
(FRBs) [4]. Radio observations of the magnetar XTE J1810—197 following its 
2018 outburst revealed a persistent 50-ms periodicity imprinted on the pulse 
profile [5]. The most commonly seen quasi-period across all observations is 
~ 76 ms (13 Hz), which is ~ P/1000. This quasi-period follows the spin-period 
scaling seen in corresponding values of the micropulses in normal pulsars [6]. 
This scaling can be most easily associated with the emission of beamlets mak- 
ing up the wider sub-pulses [7], suggesting that the periodicities are caused by 
a temporal or angular mechanism rather than the motion of the beamlets in 
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the polar cap region. Alternatively, this quasi-period could be related to sub- 
pulses or drifting sub-pulses. Each of the sub-pulses or dense, isolated ‘sparks’ 
(i.e. pair-production sites) are theorized to have a corresponding plasma col- 
umn, which radiates and generates the observed sub-pulses, which may rotate 
around the magnetic axis. Such quasi-periodic oscillations are also theorized in 
models of FRBs, where they are due to magneto-elastic axial (torsional) crustal 
eigenmodes originating close to the neutron star surface [8]. The eigenfrequen- 
cies of these oscillations are expected to depend most strongly on the neutron 
star mass and the crust equation-of-state [8]. These local crustal oscillations 
can create Alfvén waves that propagate to larger heights in the magnetosphere, 
thereby producing an oscillating F in the charge starved region to produce the 
observed coherent radio emission [9]. 

Ultimately, it is unclear what causes the quasi-periodicity in 
PSR J0901—4046. Global magnetoelastic axial (torsional) oscillations are a 
tempting explanation, but the persistence of our periodicities would require 
repeated triggers and/or very long damping times. The observed periodicities 
and frequencies, however, may be consistent with models proposed for magne- 
tars, and the similarity with the periodic feature of the radio-loud magnetar 
XTE J1810—197 are intriguing. We note that PSR J0901—4046’s position in 
the P— P parameter space is offset from the known magnetar population. We 
also note that PSR J0901—4046 may differ in other physical quantities (such 
as in its mass) which we cannot access from our observations but which are 
likely to play a role in the seismic properties of neutron stars. Hence, differ- 
ences in the behaviour compared to other neutron stars or magnetars may not 
be unexpected. It has been proposed [10] that bright coherent radio bursts 
can be produced by highly magnetized neutron stars that have attained long 
rotation periods (few 10s to a few 1000s of seconds) called Ultra-Long Period 
Magnetars (ULPMs). Recently, a source GLEAM-X J162759.5 — 523504.3 
with a period of ~ 20 minutes in the radio has been discovered, and is specu- 
lated to be a member of this class [11]. X-ray isolated neutron stars (XINS) 
are nearby cooling neutron stars with spin periods in the range 3.4-11.3 s [12] 
and are characterized by thermal, soft X-ray, emission. They are believed to 
be old, strongly magnetized neutron stars despite their non-detection in the 
radio so far [13]. A few XINS lie above the low-twist death line in Figure 1, 
implying possible ULPM origins. Interestingly, PSR J0901—4046 also falls in 
the parameter space (see Figure 1) where these ULPMs are expected to exist. 
PSR J0901—4046 could potentially be an old magnetar or a member of the 
ultra-long period magnetars, a result that needs to be confirmed with future 
multi-wavelength observations. PSR J0901—4046 is therefore an important 
piece in the puzzle of the evolution of highly magnetized neutron stars and 
their connection to FRBs. 

Typically, when magnetars are radio active, there is also often X-ray emis- 
sion. We therefore observed PSR J0901—4046 in the X-rays using Swift/XRT 
simultaneously with the MeerKAT observations on 2021-02-01 and 2021-02-02 
and did not detect any X-ray emission. Assuming a blackbody spectrum with 
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temperature 1.5 keV, and an equivalent column density of Ny = 4.32x10?! 
cm, we place 3-sigma upper limits of Lx (2 —10keV) < 1.6x10°° erg s7! 
and Ly2(2—10keV) < 3.2x 10° erg s~! on the X-ray luminosity for distances 
d, + 0.3 kpc and dz = 0.5 kpc, respectively. The location of PSR J0901—4046 
in the P — P parameter space is consistent with it having spun-down from a 
magnetar-like period of 10 s in ~ 5 Myr, assuming a braking index of 3. How- 
ever, we do not find any evidence for radical changes in the P as seen in most 
magnetars in the 7.4 months since discovery. Additionally, while magnetars 
are observed to have shallow radio spectra (e.g. [14, 15]), PSR J0901—4046 
has a measured L-band in-band spectral index of —1.7 + 0.9, which is more 
consistent with the pulsar population. Canonical, rotation powered pulsars 
are observed to have X-ray luminosities much smaller than their spin-down 
luminosities, with on average Lx ~ 1073 [16]. Conversely, magnetars are 
seen to have Ly = E. For PSR J0901—4046 , based on the X-ray luminosity 
upper limit and the spin-down F in Table 1, we see Lx Ss 102E. This places 
it closer to magnetars but is not constraining. Additionally, the single pulse 
brightness is seen to vary significantly in the 8,726 2-second integration time 
images, across the 6 L-band and 1 UHF-band epochs. The source appears to 
have secularly grown fainter (see Extended Data Figure 5 and Table 2), from a 
mean pulse brightness of 16.4 + 7.9 mJy beam~! for the observations centered 
on 59246.087481292554 to 12.9 + 5.2 mJy beam~! on 59343.62301600376, 
suggesting a dynamic magnetosphere transforming on timescales much faster 
than associated with the characteristic age, 7. If this is indeed part of a long- 
term dimming of the source, and not a short-term variation, then it is also 
reminiscent of radio-loud magnetars transitioning into quiescence. 

The single pulse polarization profiles of PSR J0901—4046 show complex 
structure, and on average are more circularly than linearly polarized (see 
Extended Data Figure 6). This is not unexpected in radio-loud neutron stars, 
particularly magnetars. The magnetar J1622—4960 exhibits different cate- 
gories of pulses of varying polarization fractions. One particular category 
shows a higher value of circular polarization. The Faraday rotation measure 
(RM) towards PSR J0901—4046 is measured to be —64+2 rad m~?. The RM 
of PSR J0901—4046 is consistent with the contribution from the smoothed 
Galactic foreground [17] and with the RMs of nearby pulsars. This there- 
fore precludes the presence of a significant intrinsic RM imparted at the 
source. A phase resolved histogram of the polarization position angle shows the 
characteristic S-shaped curve expected from a rotating magnetic dipole (see 
Supplementary Figure 3). This suggests that the line-of-sight passes close to 
the magnetic pole as we see the S-shaped curve even within a 1% duty cycle. 
This is consistent with our constraint on the impact parameter of 6 < 0.2°, 
using a rotating vector model fit. 
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4 Discussion 


The PSR J0901—4046 pulses classified as split-peak are the most common, 
~ 33% of all pulses across all observations, closely followed by a combination of 
the quasi-periodic and partially nulling pulses which together form ~ 34%. The 
normal and spiky pulses comprise ~ 27% and ~ 6% respectively. A comparison 
of the energies of the various pulse shape archetypes shows that despite the 
enormous variability seen in the pulse profile shapes, their energies span more 
or less the same range (see Supplementary Figure 8). For instance, we lose 
~ 40% of the energy to the dropouts/dips seen in the quasi-periodic and 
partially nulling pulses, which when accounted for by modelling the pulse 
envelope, is similar to the energy distribution of the ‘split-peak’ and possibly 
also the ‘normal’ pulses. This suggests that the pulses with dropouts/dips are 
not drastically brighter than the other types, implying that an overall increase 
in particle flow cannot be responsible. 

The measured period implies an extremely large pulsar light-cylinder (of 
radius Ric = cP/2n = 3.62 x 10°km) and consequently a relatively compact 
polar cap (with radius Rp = \/27R3/cP = 16.62m, where R = 10km). For an 
assumed emission height of hundreds of kilometers above the surface, the beam 
width, and consequently the duty cycle of PSR J0901—4046 is small (~1%) 
and found to be consistent with the empirical scaling relation between pulse 
width and spin period (W « P~°°) observed in canonical pulsars (e.g. [18]). 
PSR J0901—4046 is seen to lie in the P— P pulsar parameter phase space 
far from the other recently discovered slowest spinning radio emitting pulsars 
with periods of 23.5s [19] and 12.1s [20] respectively. It is also located beyond 
the ‘death line’ as defined by the RS75 [21] and CR93 [22] inner vacuum- 
gap (IVG) curvature radiation models for radio emission (see Figure 1). These 
models suggest pulsars in this region cannot support the pair-cascade pro- 
duction just above the pulsar polar cap in their inner magnetospheres that is 
required to sustain the observed radio emission. This is because at large spin 
periods, it is no longer possible to achieve the increase in thickness of the vac- 
uum gap above the neutron-star polar cap needed to maintain the required 
potential difference for pair production. This leads to the cessation of radio 
emission. However, PSR J0901—4046 does lie above the space-charge-limited 
flow (SCLF) radio emission model death line where pair cascade can be sup- 
ported through non-relativistic charges flowing freely from the polar cap if 
there is a multipolar magnetic field configuration. Unambiguous signatures of 
the presence of multipolar components close to the neutron star surface have 
been seen in magnetars [SGR 0418 +5729; 23], and more recently in an accret- 
ing millisecond pulsar [PSR J0030 + 0045; 24, 25] suggesting a likely ubiquity 
of a multipolar magnetic field configuration in neutron stars. 

The putative boundary for radio quiescence in Figure 1 indicated by Bo; 
lies about an order of magnitude below the position of PSR J0901—4046. The 
quantum process of single photon pair production (y 4 ete) is expected to 
dominate below the B,, line resulting in predominantly ‘radio loud’ pulsars. 
The quantum process of photon splitting (y > yy) is expected to dominate 
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above the B,, line resulting largely in ‘radio quiet’ pulsars due to the suppres- 
sion of pair creation. PSR J0901—4046 lies above, and at a similar distance 
from the B,, line as many of the magnetars. Unlike magnetars the radio emis- 
sion of PSR J0901—4046 has a small duty cycle, but like the magnetars it is 
highly variable. High-B radio pulsars have on occasion been observed to exhibit 
magnetar-like activity and have been termed ‘quiescent magnetars’ [26]. Radio 
emission from magnetars is usually transient and often follows a high-energy 
outburst (e.g. [27]). It is therefore useful to see how long this source has been a 
radio emitter and if any previous unidentified high-energy transient has been 
seen in this region. We did not find any historical high-energy transients coinci- 
dent with the location of PSR J0901—4046. Unfortunately, none of the relevant 
radio continuum surveys, TIFR GMRT Sky Survey (TGSS) [28], Sydney Uni- 
versity Molonglo Sky Survey (SUMSS) [29] or the Rapid ASKAP Continuum 
Survey (RACS) [30], were sensitive enough to detect the source given its cur- 
rent time-averaged flux of a couple of hundred micro-Janskys. Analyses of the 
Parkes Multibeam Pulsar Survey (PMPS) data [31] from nearby pointings also 
did not detect the source. The nearest pointing should have been sensitive 
enough, but a combination of radio frequency interference and the hardware 
high-pass filter likely prevented a detection. 

The discovery of a ~~ 117s and 118 s periodicities in the multi-wavelength 
(including radio) brightness changes of AR Scorpii (AR Sco; a radio pulsat- 
ing white dwarf binary system) [32] was interpreted as dipole emission from a 
spinning down of a magnetic white dwarf and not as a neutron star [33, 34]. 
Given the similarity in period to PSR J0901—4046 we therefore searched for 
multi-wavelength counterparts in archival data to determine whether it could 
be a related system. We identified a 17th mag Gaia source, offset by ~ 1” 
in right ascension and ~ 3” in declination from the radio coordinates, as 
a possible optical counterpart. Initial follow-up photometry with the South 
African Astronomical Observatory (SAAO) 1-m telescope showed indications 
of long term variability in the star. However, spectroscopic observations with 
the Southern African Large Telescope (SALT) revealed the optical source to 
be an A-type star, with narrow Balmer absorption lines. As we see no evidence 
for hydrogen or helium emission lines and no distinct secondary star compo- 
nent in the spectrum, we rule out the possibility of PSR JO901—4046 being an 
AR Sco type system, or associated with this A-type star. There are no other 
obvious counterparts brighter than 20 — 21 mag in this region. While the spin- 
period of PSR J0901—4046 might be consistent with a white dwarf, we do not 
see any multi-wavelength support for this. 

To ascertain if there is an un-pulsed radio component which might be 
attributed to a pulsar wind nebula, or perhaps indicate emission of a non- 
neutron star origin, we imaged follow-up MeerKAT visibility data that was 
recorded at a higher time resolution of 2 seconds. After removing the epochs 
that contain pulsed emission we obtain a 30 upper limit on the peak brightness 
of persistent radio emission of 18 wJy beam! at 1284 MHz (see Extended 
Data Figure 7). We also find no evidence for pulsed or continuum emission 
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outside of the narrow pulse window, but the radio shell is still present. The 
extreme ratio of the peak on-pulse flux to the off-pulse flux, the large first 
period derivative, the timing properties, and the lack of evidence for detections 
at other wavelengths supports our hypothesis that PSR J0901—4046 is a radio 
emitting neutron star with one of the longest known periods. 


5 Implications for the population of 
radio-emitting neutron stars 


Although modern pulsar surveys are sensitive to a wide range of radio emitting 
neutron stars, the serendipitous discovery of PSR J0901—4046 has revealed 
some of the biases that still remain, and highlights that there may be many 
more sources like this to be found. The long duration of the pulse, low DM and 
long period are problematic for commonly used single pulse and periodic pulsar 
search techniques. The very narrow duty cycle suggests a strong bias where 
many other similar systems may have beams missing the Earth completely. 
This suggest that there are many more neutron stars in the Galaxy than the 
known population suggests, unless many pulsars continue to emit for longer 
than previously thought or if there is an evolutionary link to another class of 
neutron star such as the magnetars [35], or perhaps a combination of all of 
these. The position of PSR J0901—4046 in the P — P parameter space along 
with the unusual single pulse properties such as quasi-periodicity and partial- 
nulling, make it a potentially very useful target for understanding the radio 
emission properties of neutron stars across the population. Future image and 
time domain searches for similar long-period objects could prove vital to our 
understanding of the Galactic neutron star population and potentially links to 
FRBs. 
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Calibration of interferometric imaging data 


The MeerKAT observations of the field around PSR J0901—4046 are sum- 
marised in Table 2. Nine distinct observations were used, eight of which 
used MeerKAT’s L-band (856 — 1712 MHz) receivers, and one of which was 
observed at UHF band (580 — 1015 MHz). Total on-target times (Tops) and 
correlator integration times per visibility point (Tint) are listed in Table 2. 
The latter is what limits the time resolution of any individual image of the 
field. The ‘discovery’ observations, associated with observations of Vela X-1 
from the ThunderKAT project, were taken with the correlator configured to 
deliver 32,768 spectral channels. The follow-up (DDT) observations that tar- 
geted PSR J0901—4046 directly used 4,096 channels. For the imaging data the 
observations were in all cases averaged down to 1,024 channels prior to the 
commencement of the processing. 

The approach to imaging the field was common for all observations. Each 
observation contains 5 minutes scans of the standard primary calibrator source 
J0408—6465, and the scans of the target field were bracketed by observations of 
the nearby secondary calibrator J0825—5010, which was observed for 2 minutes 
for every 15 minutes on the target for the ThunderKAT data, and for every 30 
minutes for the DDT observations. The calibrator scans were flagged in order to 
remove radio frequency interference, and the low-gain edges of the telescope’s 
bandpass response. Bandpass, delay, and flux-scale corrections were derived 
from the observations of the primary calibrator, and time-dependent complex 
gain and delay corrections were derived from the scans of the secondary. These 
corrections were then applied to the target data. These steps were all performed 
using the CASA package [36]. 

Following the application of the referenced calibration, the target data were 
flagged using the TRICOLOUR (https: //github.com/ska-sa/tricolour/) software. 
The target field was imaged using WSCLEAN [37]. Deconvolution was allowed 
to proceed in an unconstrained fashion. The field exhibits some complex radio 
morphology, thus a cleaning mask was derived from the first image, after 
which the imaging was repeated with deconvolution only proceeding within 
the masked regions. The frequency dependence of the sky was captured by 
imaging the data in eight separate sub-bands, using a fourth-order polynomial 
fit to capture spectral curvature, mainly an instrumentally-induced property 
due to the frequency dependent antenna primary beam response and the broad 
bandwidth. A sky area (3.12 x 3.12 deg?) much larger than the main lobe 
of the primary beam (~1 deg at FWHM) was imaged in order to decon- 
volve bright off-axis sources that are detected through the primary beam 
sidelobes. The use of the cleaning mask, spectral settings, and large sky area 
in this second imaging run are all to ensure a reliable model for subsequent 
self-calibration, which consisted of solving for instrumental phase and delay 
corrections for every 32 seconds of data using the CUBICAL package [38]. The 


Springer Nature 2021 TeX template 


A radio-loud neutron star with an ultra-long spin period 11 


scripts used to perform the data reduction process also provide an exhaus- 
tive list of the calibration and imaging parameters, and can be found online 
(https://github.com/IanHeywood/oxkat v0.2) [39]. 


Snapshot imaging of PSR J0901—4046 


To expedite the production of per-integration time images, we first subtract a 
model of the sky that captures most of the bright emission, but critically does 
not include any clean components that are associated with PSR J0901—4046 
itself. For each MeerKAT observation, the self-calibrated data were imaged, 
the resulting model images were masked at the position of PSR J0901—4046 , 
and the modified model was inverted into a set of model visibilities, which were 
then subtracted from the data. Images could then be made for each correlator 
dump time (8 or 2 seconds) to search for pulsed emission from the target. In 
the case of the ThunderKAT data the visibilities were first phase-rotated to the 
position of PSR J0901—4046. Since the dominant emission in the field has been 
subtracted, small images around the target are viable, and no deconvolution 
needs to be performed under the (valid) assumption that PSR J0901—4046 is 
unresolved by MeerKAT. This speeds up the imaging process considerably. 

Extended Data Figure 5 shows the peak brightnesses of the pulses as 
detected in 8,726 2-second snapshot images from the six L-band epochs with 
this integration time. The pulse brightness varies significantly, however no 
pulses are missed at the sensitivity limit of our observations, with the excep- 
tion of the cyan region in the lower left panel, where data were lost due to 
lightning. The mean RMS noise the snapshot images being 350 wJy beam™! 
with a standard deviation of 50 uJy beam~!. The right hand column of panels 
in Extended Data Figure 5 shows the pulse brightness expressed as a signal- 
to-noise (S/N) ratio. The brightness is measured as the peak pixel value in a 
400 pixel box centered at the position of PSR JO901—4046. The noise is taken 
to be the standard deviation of the pixels in an off-source box of equivalent 
size. The blue curve on the right hand column of panels shows the S/N ratio 
of the peak pixel in the off-axis box. 


A search for persistent (off-pulse) radio emission 


Jointly imaging and deconvolving the visibilities used to produce the 2-second 
images results in the image shown in the left panel of Extended Data Figure 7. 
The accumulated pulse emission results in the prominent compact source in the 
centre of the image, with a peak brightness of 40 (45.2) wJy beam~!. Identify- 
ing the timestamps of the pulses shown in Extended Data Figure 5 allows us to 
re-image the data with those integration times excluded in order to search for 
off-pulse (persistent) radio emission associated with PSR J0901—4046. This 
process results in the image shown on the right hand panel of Extended Data 
Figure 7. There is a 4.3 pJy beam~! (~10; ¢ = 4.7 Jy beam~!) peak in the 
pulse-subtracted radio map spatially coincident with the peak of the pulsed 
emission in the image formed from the full dataset. We can thus place a 30 
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upper limit on the peak brightness of a persistent radio source coincident with 
the peak of the pulsed emission of 18 pJy beam~!. The other sources visible in 
Extended Data Figure 7 are faint compact sources that were not deconvolved 
in the per-epoch imaging process, and thus not subtracted from the visibilities. 

In the deepest image we have made in Extended Data Figure 7, there is 
diffuse emission in the region of PSR J0901—4046. More analysis is needed to 
see if it is somehow associated with PSR J0901—4046 as this is a complex 
region of the sky with lots of diffuse emission. However if it were attributable 
to radio emission from a supernova remnant that could be associated with 
PSR J0901—4046 it would suggest that the source was much younger than the 
characteristic age and have important implications for its evolution. None of 
the features visible in Extended Data Figure 7 are above the noise floor of the 
2 second images, and thus do not contaminate the measurements presented in 
Extended Data Figure 5. 


DM estimate 


The DM of each single pulse was estimated by maximizing for structure within 
the burst envelope using DM_phase (https://github.com/danielemichilli/DM_ 
phase). The structure-optimized DM is determined by maximizing the coher- 
ent power across the bandwidth [40]. We de-dispersed the data over a trial 
DM range of 49.0 < DM < 54.0 pe cm~® in steps of 0.1 pe cm~?. The 
uncertainty on each DM estimate was calculated by converting the standard 
deviation of the coherent power spectrum into a standard deviation in DM via 
the Taylor series. We measure a weighted average DM of 52+1 pccm~? for 
PSR J0901—4046. 


Timing 

The MeerTRAP pipeline searches data in real time and for each transient 
event detected, it writes out a short SIGPROC filterbank file that contains a few 
seconds of the original input data stream centered around the detection time of 
the associated event. For each detection, a substantially smaller, second-stage 
candidate file is also created as follows: the native resolution filterbank file is 
dedispersed at the detection DM reported by the search pipeline, a reduced 
time span of the data window equal to the dispersion delay of the candidate 
DM is extracted and lastly the time and frequency resolution of the data 
are reduced to an appropriate level according to the reported pulse width for 
the event (larger widths correspond to a larger acceptable degradation factor 
of the data). Second-stage candidate files are small enough to be stored en 
masse, but the native resolution filterbank files are not; only those deemed 
very likely to contain a genuine astrophysical event by an automated classifier 
are kept, the rest are otherwise regularly deleted. For the original detection 
of PSR J0901—4046, we thus had access to a filterbank file with our native 
resolution of approximately 306.24 ys and 1024 channels across the 856 MHz 
bandwidth at L-band. A detailed visual inspection of second-stage candidate 
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plots around the time of discovery showed that we had actually made more 
detections of the source, but the unusual nature of the wide pulses and the 
relatively low DM meant that they were not initially labelled as astrophysical 
and the associated filterbank data had already been deleted. However, there 
is still sufficient information in the second-stage candidate files to allow the 
manual determination of an arrival time, albeit with significantly increased 
uncertainty; we thus obtained a total of 14 approximate arrival times from the 
day of the discovery. 

The uncertainties on those times of arrival (ToAs) were estimated as part 
of the initial periodicity analysis, which consisted of a simple Bayesian linear 
regression with three parameters: a pulse period, an initial phase term, and 
an uncertainty scaling factor where the underlying assumption was that the 
uncertainty on each arrival time was proportional to the pulse width reported 
by the search pipeline. This yielded the initial period estimate of P = 75.89+ 
0.01 seconds previously mentioned, and a root mean square uncertainty of 100 
ms on the arrival times. These 14 initial ToAs were then used along with the 
imaging data to further constrain the pulse period. They are also included in 
the timing analysis (see the orange points in Extended Data Figure 1). 

We generated ToAs for each of the two 30 minutes observations on each 
of the 6 MeerKAT L-band observing epochs between February and May 2021. 
The filterbank data recorded with the TUSE instrument were used for tim- 
ing. These data also had a resolution of approximately 306.24us and 1024 
channels across the 856 MHz bandwidth at L-band. The data were folded on 
to 65536 phase bins and incoherently dedispersed using the DSPSR software 
package [41]. Further manipulation of the data used the tools available in the 
PSRCHIVE package [42]. Radio-frequency interference was removed manually 
using PSRZAP. 

To maintain a coherent timing solution, we secured available observing slots 
at the Parkes radio telescope in Australia and monitored PSR J0901—4046 
over 4 epochs. The observations were performed using the Ultra-Wideband 
receiver (UWL) across a 704-4032 MHz band. The data were coherently 
de-dispersed at the DM given in Table 1, divided into 1024 phase bins for 
each of 3328 frequency channels of 1 MHz bandwidth, and written to disk. 
The duration of the observations varied between 1 and 2 hours. The RFI 
environment at Parkes is such that many frequency channels in the data 
contained strong baseline fluctuations on a time scale comparable to the 
typical pulse duration of PSR J0901—4046, which required applying addi- 
tional layers of RFI mitigation to make the averaged pulse unambiguously 
identifiable in each integrated observation. An updated version of the clfd 
(https://github.com/v-morello/clfd) RFI cleaning package described in [43] 
was used to that effect. The relatively steep spectrum of the source meant 
that it was not easily detectable above 1.8 GHz even after cleaning and so we 
only used frequencies in the range 0.704 to 1.8 GHz GHz for generating arrival 
times. 
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The initial folding and timing analysis for both telescopes used the best 
known period and DM at the time and the position determined from the 
imaging. For the MeerKAT data a noise-free template was made by fitting 
von Mises functions to the data from the first epoch in February 2021 using 
the PSRCHIVE program PAAS. A similar procedure was followed for the Parkes 
UWL data also using the first observation from February 2021. For both the 
MeerKAT and Parkes data ToAs were obtained by using PSRCHIVE’s PAT 
which cross-correlated the template with an average profile for each of the 30 
minute observations. A ToA (the first red diamond in Extended Data Figure 
1) was also determined from the single pulse obtained from the filterbank file 
saved from the discovery observation by cross-correlating it with the MeerKAT 
template (As the filterbank file was significantly shorter in duration than the 
period of the source we padded it with random Gaussian noise so that it could 
be folded using DSPSR and preserve the timing information). 

Timing was done using TEMPO2 [44] with the JPL DE436 planetary 
ephemeris _— (https: //naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de435s. 
bsp.Ibl). The ToAs were fit using a model including the period P and period 
derivative P and a jump between the UWL data and the MeerKAT data 
(Another, fixed jump of 6.115060037383178s, was needed for a few of the 
arrival times from the discovery epoch due to an uncertainty of one data 
buffer that occurred in these early data). We do not fit for position as it is 
well determined from the imaging as described in previous sections. An astro- 
metric precision of about ~ 1” results in a change in arrival time over a year 
of less than 2.3 ms, so is not affecting the measured parameters, especially 
the P which can show a covariance when there is less than a year of data. We 
also do not fit for DM as this is sufficiently well determined from optimizing 
the S/N of the individual pulses and a jump is fit between the UWL and 
MeerKAT data. 

Once a coherent timing solution was obtained across the entire data set 
the filterbanks were refolded and dedispersed and a new noise-free template 
was made, based on the sum of all of the detected pulses, and new ToAs were 
obtained and a new final timing solution was determined and is presented in 
Table 1 with the residuals shown in Extended Data Figure 1. Each 30-minute 
observation is the average of about 24 pulses and so there is some pulse phase 
jitter, which can be seen in the MeerKAT data, where the error bars are 
smaller than the data points and also less than the scatter in the arrival times. 
However the overall timing RMS is only 5ms which is just under 1/10000th 
of the pulse period which is approaching that of the best millisecond pulsars 
and is attributable to the very high signal to noise but also suggests that the 
arrival times are not significantly affected by the pronounced variations in the 
pulse properties and likely reflects the overall similarity of the pulse envelope. 
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Quasi-periodicity Analysis 


Radio loud neutron stars are seen to exhibit a rich variety of intensity variations 
over timescales of microseconds to years. Within an individual pulse, sub- 
structures manifest most conspicuously as sub-pulses/components that have 
random-like but also, pulse- and source-dependent continuous or quasi-periodic 
variability in time intervals. Quasi-periodicities are usually seen as repeat- 
ing “micropulses” [45] forming part of “microstructure” superimposed on the 
wider sub-pulses (e.g. [46]). Microstructure is usually theorized to be caused by 
mechanisms related to magnetospheric radio emission. There is often a variety 
of timescales observed, even within a given source (e.g. [47]), for pulsars with 
typical periods of about ~ 1s, sub-pulses tend to have widths of a few to tens 
of ms. Timescales and periodicities in the shorter micropulses tend to scale 
like P,, ~ P/1000 [6, 7], where often the micropulse duration and periodicity 
scale like w, ~ P,,/2 [6]. 

We followed standard methods [6, 47] to determine the timescales of the 
short-time structure for PSR J0901—4046. An auto-correlation function (ACF) 
of pulses detected by the TUSE pipeline containing both sub-pulses and 
microstructure, will show a peak at zero-lag corresponding to the DC compo- 
nent, followed by a peak at short timescales due to possible microstructure and 
then a second peak associated with the sub-pulse structure [48]. We extracted 
the timescale of the sub-pulses by performing an ACF analysis of the single 
pulse intensities. We compute the cross-correlation of the de-dispersed signal 
as a function of time with a delayed copy of itself given by, 


acr(r) = f fH) f(t — 7) dt, (1) 


where 7 is the time lag. The zero-lag value, associated with self-noise, was 
excised from the ACF. Given the complexity of the structures seen in the single 
pules, we visually inspected each pulse to understand what can and cannot be 
inferred about the timescales. The characteristic separation or quasi-period of 
the sub-pulses (defined as P)) measured across the whole observing band is 
given by the time lag of the peak of the first feature following the zero-lag in 
each ACF. We only measure the the quasi-periods for those pulses which are 
visually obvious in the ACFs. 

As a result, we observe some of the quasi-period values to be harmon- 
ically related as shown in Extended Data Figures 2 and 4, and in some 
cases, the separations between the peaks is almost the separation between the 
dips or dropouts in power. Occasionally, two or more quasi-periods co-exist 
in a single rotation. Upon visual inspection we notice that some pulses in 
PSR J0901—4046 exhibit variations in widths as well as quasi-periods within 
a single rotation as seen in Extended Data Figure 3. We do not observe these 
quasi-periodic pulses to follow a trend in time, nor do they appear to precede 
or follow any other particular type of pulse. 
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These quasi-periodic features could be interpreted as sub-pulses or drifting 
sub-pulses. Although the latter are usually characterized by a fixed, from pulse 
to pulse, separation between sub-pulses. The sparking discharge from IVG 
models explains sub-pulses through non-stationary spark-associated plasma 
flows. The IVG is discharged in the form of dense isolated sparks (i.e. pair- 
production sites), where the lateral size and distance between the sparks are 
comparable to the gap height [21]. Each spark has a corresponding coherent 
plasma column, which radiates and generates the observed sub-pulse compo- 
nents in pulsars. According to this ‘Spark’ model [49], up to three sparks can 
be accommodated in the polar cap of PSR J0901—4046 , which is much fewer 
than the number of features typically seen in its single pulses. 

The dominant quasi-period of ~ 76 ms (corresponding to a frequency of 
13 Hz) follows the spin-period scaling seen in corresponding values of the 
micropulses in normal pulsars [6] (see Extended Data Figure 4). This scal- 
ing can be most easily associated with the emission of beamlets making up 
the wider sub-pulses [7], implying that the periodicities are caused by a tem- 
poral or angular mechanism rather than a radial mechanism. Apart from the 
consistent scaling in the value of the quasi-period, the appearance of sev- 
eral periodicities within one sub-pulse has also been seen in normal pulsars 
(e.g. [6, 45, 47]). Overall, this makes it tempting to associate these structures 
to “normal” microstructure in pulsars in a similar manner. 

However, the appearance of the dropouts (see e.g. the top-left example in 
Extended Data Figure 2) is different to that of normal micropulses. In con- 
trast, it is very reminiscent of quasi-periodic oscillation (QPO) features seen 
both in the emission of hard short X-ray bursts and the tail of energetic giant 
flares of magnetars. The “dropout” pattern seen in the quasi-periodic and par- 
tially nulling pulses is very unusual for pulsar radio emission. Nevertheless, 
we establish that these dropouts are a genuine feature of the emission of the 
source. We see these features at all epochs in the filterbank data recorded 
by both the TUSE and APSUE instruments, as well as the raw voltage data 
extracted directly from the F-engine. These features track the dispersion of 
the source and, importantly, are not seen in the “off pulse” emission. Addition- 
ally, we do not observe these features in any of the test pulsar (JO820—4114) 
observations at the start of each epoch, thereby ruling out the possibility of 
instrumental artifacts. The magnetar QPOs with frequencies between 18 to 
1800 Hz are often interpreted as a result of seismic vibrations of the neutron 
star (see e.g. [26, 50] and references therein). Global magnetoelastic axial (tor- 
sional) oscillations are expected to be able to explain the frequencies as low as 
seen here. They have also been put forward as explanations for sequences of 
emission features seen in Fast Radio Bursts (FRBs) invoking magnetar oscilla- 
tions [8]. Adopting such an explanation for the quasi-periodicities observed in 
PSR J0901—4046 would demand, given their persistence in our observations, 
a repeating triggering mechanism or modes with long-lived eigenfrequencies. 
The nature of the modes, their eigenmodes and, importantly, their damping 
times depend strongly on the physics and properties of the neutron star’s crust, 
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the mass, the equation-of-state and, to some degree, also on magnetic field 
strength [26]. Most references discuss damping time lengths only for the dura- 
tion of the seen QPOs or FRB emission sequences, i.e. tens of ms to seconds 
(e.g. [8]), which is clearly too short to explain the consistency of our observed 
periodicities over many weeks and months (see Extended Data Figure 4). 
Interestingly, radio observations of the magnetar XTE J1810-197 during its 
renewed radio brightness following its 2018 outburst also revealed a persistent 
50-ms periodicity in its pulse profile seen for about 10 days. As reported by 
Levin et al. (2019) [5], this emission feature was imprinted on the pulse profile 
simultaneously at a range of observed radio frequencies between 1 and 9 GHz. 
The frequency of this feature of about 20 Hz is obviously similar to that of 
our dropouts in PSR J0901—4046. In XTE J1810—197 the feature showed a 
remarkable constancy in phase relative to the main pulse profile, implying that 
the periodicity is not a temporal modulation of the emitting source, but must 
be due to a periodic structure in the radiation beam pattern that sweeps across 
the Earth as the pulsar rotates. Levin et al. suggested that this pattern could 
arise from a stable structure on the surface of the neutron star at the base of 
the magnetic field lines hosting the emitting particles for the radio component. 
The stability of the pattern would require a frozen-in wave pattern of radial 
dimension of the surface height, temperature or magnetic field. Such a pattern 
would be reminiscent of surface waves in the neutron star crust, similar to 
those discussed above, but they would need to be stable over at least 10 days. 


X-ray Follow-up 


We requested Neil Gehrels Swift Observatory (Swift) observations to search for 
a candidate X-ray counterpart to PSR J0901—4046. We obtained three obser- 
vations for a total exposure of 7419.030 s (ObsIDs 00014019002, 00014019003, 
00014019004, taken on MJD 59245.83, 59352.04, and 59353.29 with individ- 
ual exposures of 3872.021 s, 1345.595 s, and 2703.971 s, respectively). The 
three observations were separated by 1, 9, and 10 days from their closest 
radio observation, respectively. We extracted an image using the XRT product 
generator online reduction pipeline (https://www.swift.ac.uk/user_objects/ 
index.php) [51]. A visual inspection of the image showed that the source 
was not detected, therefore we use the SOSTA tool within the XIMAGE 
environment to perform source statistics. SOSTA allows one to use a local 
background to determine the significance of a source and its count rate, 
rather than a global background estimate. Extracting events at the nominal 
position of PSR J0901—4046 from a square box with 16 pixels side (+38’) 
we obtained a 3-sigma upper limit to the count rate of RY? = 1.57x1073 
cts/s in the 0.5—10 keV energy range. We then assumed a blackbody spec- 
trum with temperature 1.5 keV, and an equivalent column density of Ny = 
4.32107! cm~3 (i.e. the Galactic equivalent column density in the direction of 
PSR J0901—4046). Using WebPPIMS (https: //heasarc.gsfc.nasa.gov/cgi-bin/ 
Tools/w3pimms/w3pimms.pl, PIMMS v4.11b.) we estimated a 3-sigma upper 
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limit to the 0.5—-10 keV X-ray flux of F < 1.2x10713 erg cm~? s~!, which at 
a distance of d, * 328 pc and dg & 467 pc corresponds to a 3-sigma upper 
limit to the X-ray luminosity Lx, < 1.6x10°° erg s~! and Lxg < 3.2x10*° erg 
s~', respectively. 

We note that the XMM-Newton archive includes two archival observations 
of the Vela X-1 field (ObsID 0406430201 and 0841890201, taken in 2018 and 
2019, respectively). In such pointings PSR J0901—4046 is located at the very 
edge of the EPIC-MOS image, and in one of the two observations only one of 
the MOS cameras was active (the other was switched off for telemetry reasons). 
Given that the response of the instrument is not ideal so close to the edge of 
the CCD, and calibration might not be reliable, we decided to not use these 
observations, and rely solely on the more conservative - but likely more robust 
- upper limit derived from the Swift data. 
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Table 1: Pulsar timing and model parameters for PSR J0901—4046. This 
includes the measured quantities and the derived quantities from the timing 
analysis over the span of this observing campaign. Uncertainties in parentheses 
as l-o errors on the last significant quoted digit. 


Data and model fit quality 
Modified Julian Date (MJD) range............. 59119.0 to 59343.6 (7.4 months) 
Number of TOAS: a 4cccs csv ewne eee tin sate se veres 29 
Weighted root mean square timing residual (ms) 5.7 
Measured quantities 


Right ascension, a (J2000)..................008 09”017™ 29.249° + 1.0” 
Declination, 5 (J2000).................0..02000- —40°46'02.984” + 1.0” 
Pulse frequency, Vo... 6.0. cece cence cece e ences 0.013177739873 + 9.9 x 10-12 s—! 
First derivative of pulse frequency, v ........... —3.9+0.2 s~? 
Pulse period; Press aie tod iaoodaienn oe bi pewter cs 75.88554711 + (6 x 10-8) s 
Period derivative, P ........... 000 c cee eeeee eee (2.25 + 0.1) x 10-13 s s~t 
Dispersion measure, DM ..................00005 52+ 1pecm=3 
Full width at half maximum, Ws59 (L-band) ..... 299+ 1 ms 
Full width at half maximum, Ws0 (UHF-band) 296 + 4 ms 
Spectral index, @ 10... se eee eee eee eee —1.7+0.9 
Rotation measure, (RM) .................0.0008 —64+4+2 rad m~? 
Fractional linear polarization .................. 12.2+0.2 % 
Fractional circular polarization ................ 21.04+1.9% 

Inferred quantities 
Distance (YMW16), dy ........ 0.0... eee eee eee 328 pc 
Distance (NE2001), dg ............ ce cece eee 467 pc 
Characteristic age, T 00... 6... cece eee 5.3 Myr 
Surface dipole magnetic field strength, B ...... 1.3x104G 
Spin-down luminosity, BE ....... 60.0... cece eee 2.0 x 1078 erg s~! 
Period-averaged radio luminosity, Li4o0 at dz... 89 Jy kpc? 


X-ray Luminosity, Lx (0.5 — 10 keV) at dg ...... <S 3.2 108° erg s~! 
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Table 2: MeerKAT observations of the PSR J0901—4046 field. The first three 
rows labelled TKAT are discovery observations targeting the Vela X-1 field, 
while the rest labelled DDT are follow-up observations. See the text for details. 


Date Block ID RA Dec Band Nant Tabs Tint Origin 

UT, J2000 J2000 J2000 h s 

2020-09-25 1600995961 09"0206.86° —40°33'16.9" LL 59 0.5 8 TKAT 

2020-09-27 1601168939 09"02""06.86°  —40°33'16.9" LL 61 0.5 8 TKAT 
- ve 

2020-10-11 1602387062 09"027"06.86° —40°33'16.9° LL 60 0.5 8 TKAT 

2021-02-01 1612141271 09"01'"29.35°  —40°46'03.6 LL 64 1 2 DDT 

2021-02-02 1612227667 09"01'29.35°  —40°46'03.6 LL 61 1 2 DDT 

2021-02-10 1612994791 09"01™29.355 —40°46'03.6" LL 62 1 2 DDT 

2021-03-03 1614794470 09"01'"29.35° + —40°46'03.6 LL 63 1 2 DDT 
Ff ica 

2021-04-02 1617367872 09"01'"29.35°  —40°46'03.6 LL 63 1 2 DDT 
a 

2021-04-02 1617376889 09"0129.355 —40°46'03.6 UHF 62 1 2 DDT 
¥ wt 

2021-05-09 1620567645 09"01'"29.35°  —40°46'03.6 LL 62 1 2 DDT 
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Fig. 1: P — P diagram based on the ATNF pulsar catalog. The various sub- 
classes of pulsars are represented by the markers in the legend. The longest spin 
period radio pulsars and the white dwarf binary system AR Sco are highlighted 
in red. Lines of constant age and magnetic field are shown as dotted and dashed 
lines respectively. The lower right corner of the figure represents the ‘death 
valley’ with various death lines from the literature, where sources below these 
lines are not expected to emit in the radio. The solid death line represents 
Equation 9 in CR93 [22]. In dot-dashed and dashed are the death lines modeled 
on curvature radiation from the vacuum gap and SCLF models as shown by 
Equations 4 and 9 respectively in [52]. Sources above the low-twist death line 
are potential ultra-long period magnetars. 
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Fig. 2: Gallery of the pulse morphology types of PSR J0901—4046. The mor- 
phological type is given in the title for each panel. The top panels are pulses 
observed in the UHF-band while the bottom panels are pulses observed at 
L-Band. 
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Extended Data Figure 1: Timing residuals of PSR J0901—4046. The resid- 
uals from the best fit timing model given in Table 1. The orange data 
points are determined from the original MeerTRAP detection images, the first 
red diamond corresponds to a single pulse and the remaining red diamonds 
are determined from each of the half hour long follow-up observations with 
MeerKAT. The error bars are 1-a0. We used the L-band MeerKAT data for 


the timing analysis. The light coloured data points are from the Parkes UWL 
observations. 
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Extended Data Figure 2: Examples of quasi-periodic pulses. The top two 
rows show pulse profiles and their corresponding ACF's at 306.24,s resolution, 
respectively. The value the of quasi-period is indicated by the black vertical 
lines. The bottom two rows show the off-pulse regions and their corresponding 


ACFs. 
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Extended Data Figure 3: Example of a pulse exhibiting more than one 
quasi-period. Some quasi-periodic pulses as shown here, exhibit multiple quasi- 
periods within a single rotation. 
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Extended Data Figure 4: Estimates of the quasi-period across all epochs. 
The (orange) circles are the measured quasi-periods for each single pulse. The 
most commonly observed average quasi-period is 75.82 ms with the minimum 
period being 9.57 ms. The lags are arranged in lag length and not in time order. 
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Extended Data Figure 5: Radio light-curves of PSR J0901—4046. A regu- 
lar series of pulsed emission detected in the L-band snapshot imaging for six 
observing epochs. Please refer to Section 5 for details. 
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Extended Data Figure 6: Polarization profiles of PSR J0901—4046 at 
1.3 GHz and 700 MHz. Top Panel: Time series of two single pulses of 
PSR J0901—4046 at 1284 MHz. Bottom Panel: Two different single pulse time 
series at 737 MHz. For both panels, the total intensity is represented by the 
black solid line, the red solid line denotes the linear polarization while the blue 
solid line denotes circular polarization. The polarization position angle is not 
absolutely calibrated at 737 MHz. 
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Extended Data Figure 7: MeerKAT image of the PSR J0901—4046 region 
at 1.28 GHz. The left hand panel shows the image with the pulsed emission 
included, and the right hand panel shows the same field following the removal of 
the integration times containing pulses. No persistent radio source is associated 
with PSR J0901—4046 to a 3o limit of 18 pJy beam. The diffuse shell-like 
structure that surrounds PSR J0901—4046 is partially visible, possibly the 
supernova remnant from the event that formed the neutron star. 
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MeerKAT description and configuration for 
transient searches 


The Meer(more) Karoo Array Telescope (MeerKAT) [53-55] operated by the 
South African Radio Astronomy Observatory (SARAO) comprises 64, 13.5-m 
antennas distributed over 8-km in the Karoo region in South Africa. 40 of these 
dishes are concentrated in the inner ~ 1-km core. The MeerTRAP project is a 
commensal programme to search for pulsars and fast transients whilst piggy- 
backing on the large survey programmes at MeerKAT. The ThunderKAT 
programme is the image-plane transients search programme at MeerK AT and 
operates in synergy with MeerTRAP for blind transient searches. In the obser- 
vations presented in this work, MeerKAT operated at a centre frequency of 
1284 MHz with a usable bandwidth of ~ 770 MHz in the L-band observations 
and at a centre frequency of 798 MHz with a usable bandwidth of ~ 435 MHz 
in the UHF-band observations. MeerK AT observes simultaneously in incoher- 
ent and coherent modes using the MeerTRAP backend. In the coherent mode, 
the voltages from the inner 40 dishes of the ~ 1-km core of the array are 
coherently combined to form up to 780 beams on sky with an aggregate field- 
of-view of ~ 0.4 deg~?. In the incoherent mode the intensities of all available 
(up to the maximum 64) MeerKAT dishes are added to create a less sensitive 
but much wider field-of-view of ~1.3 deg?. 

The MeerTRAP backend is the association of two systems: the Filter- 
bank and Beamforming User Supplied Equipment (FBFUSE), a many-beam 
beamformer that was designed and developed at the Max-Planck-Institut fiir 
Radioastronomie in Bonn [56, 57], and the Transient User Supplied Equipment 
(TUSE), a real-time transient detection instrument developed by the Meer- 
TRAP team at the University of Manchester [58]. The input data stream to 
FBFUSE consists of the complex-valued channels from every dish produced by 
the MeerKAT F-engine, which are transmitted via the Central BeamFormer 
(CBF) network; FBFUSE applies the geometric and phase delays (obtained 
by observing a bright calibrator) before combining the data streams from the 
dishes into one incoherent beam and up to 780 total intensity tied-array beams. 
The beams can be placed at any desired locations within the primary beam of 
the array, but are by default tessellated into a circular tiling centered on the 
boresight position, and spaced so that the response patterns of neighbouring 
beams intersect at the 25% peak power point. The beams are then sent over 
the network to TUSE for processing. 


Transient search pipeline 


TUSE consists of 67 Lenovo servers with one head node and 66 compute nodes. 
Each compute node contains two Intel Xeon CPUs with 16 logical cores each, 
two Nvidia GeForce 1080 Ti Graphical processing units (GPUs) and 256 GB 
of DDR4 Random Access Memory (RAM). Each of the nodes is connected to 
a breakout switch via 10 GbE network interface cards (NIC) that are used to 
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ingest data from FBFUSE. The data from FBFUSE are received on the NICs 
as SPEAD2 (https: //casper.ssl.berkeley.edu/wiki/SPEAD) heaps: small, self- 
describing units that contain a header and 8KB of data in time-frequency 
order, attached to a specific frequency sub-band and beam. FBFUSE streams 
these heaps over a range of multicast network addresses, where each address 
is associated with a fixed subset of 12 coherent beams in a typical configura- 
tion. Each TUSE node then subscribes to a distinct multi-cast address, which 
effectively distributes the processing load evenly across the TUSE cluster. The 
data ingest code on a TUSE node collects heaps and writes their contents to 
POSIX shared memory ring buffers (one per beam), so that the data for each 
beam are reconstructed in their natural time-frequency order (i.e. with the 
frequency axis contiguous in memory). Once a data segment has been fully 
reconstructed, it is ingested by the search pipeline. More details on TUSE will 
be presented in an upcoming paper. 

We utilize the highly-optimized Graphics Processing Unit (GPU)-based 
ASTROACCELERATE software [59, 60] to search for dispersed signals. The real- 
time search is performed by incoherently de-dispersing in the DM range 
0-5118.4 pe cem~? at L-band and 0-2664.9 pc cm~? in the UHF-band, divided 
into multiple sub-ranges with progressively increasing DM steps and time aver- 
aging factors. We also search up to maximum boxcar widths of 0.67 s and 
0.77 s in the L-band and UHF-band respectively. The extracted candidate files 
contain raw filterbank data of the dispersed pulse and additional padding of 
0.5 s at the start and at the end of the file. 

For the targeted observations, where the position of PSR J0901—4046 is 
known to within an arcsecond, we were able to run in a mode where only 2 
nodes with one beam per node were processing the data in real time. Only 
extracted candidates were saved for further examination. In addition to the 
real-time processing, two more nodes were used to record the continuous data 
to disk at a full data rate for further offline processing. 


APSUSE data recording 


The APSUSE instrument consists of 68 high-performance Huawei FusionServer 
2288H v5 servers with two head nodes for control. Nodes are split into two 
classes depending on their primary use. There are eight ingest nodes, designed 
for the capture and collation of Ethernet data and its writing to disk. Each 
ingest node has two 3-GHz Intel Xeon Gold 6136 CPUs, a 40-GbE Ethernet 
adapter and an 56 Gb/s FDR Infiniband adapter. The remaining 60 nodes 
are used for data analysis, with each having two 2.1-GHz Intel Xeon Silver 
4116 CPU and two Nvidia GTX 1080 Ti GPUs. Additionally these nodes 
have FDR Infiniband adapters. All nodes in the cluster are connected via an 
Infiniband fabric with each node hosting 8 x 8 TB hard drives, providing a 3.2- 
PB distributed file-system with 40 GB/s of write performance. The APSUSE 
ingest nodes were used during PSR J0901—4046 observations to record beams 
formed by the FBFUSE instrument. At L-band, data were recorded with 4096 
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frequency channels over a 856-MHz band centered at 1.284 GHz with a time 
resolution of 76.56 ys. At UHF the data were recorded with 4096 frequency 
channels over a 544-MHz band centered at 816 MHz with a time resolution of 
120.47 ws. 


Full-time resolution data capture pipeline 


We also used the transient buffer mode of the MeerTRAP pipeline to capture 
complex channelized data in order to investigate the sub-pulse structure and 
polarization properties of the PSR J0901—4046 single pulses. For each obser- 
vation, the MJD of the first detection in the real-time pipeline was used as a 
reference in order to predict the time of arrival of future pulses using the best 
known pulse period. Then, a list of predicted MJDs were sent as a trigger to the 
transient buffer capture code (https://github.com/ewanbarr/psrdada_cpp) to 
each beamformer node. On receiving the trigger, the code uses the DM of the 
source to extract the complex channels around the time of the pulse account- 
ing for the dispersion delay in each frequency channel before saving the data 
to disk. Each beamformer node saved data for one sub-band and all antennas 
used in the observation. For each pulse, around 1.5 seconds of complex channel 
data were stored. For every pulse, the data for each sub-band were combined 
and the proper complex gain correction applied per antenna and per frequency 
channel. The complex channels from each antenna were then added coherently, 
in phase, to form a beam at the phase centre of the telescope. These coherently 
beamformed complex channels were then used to study the polarization prop- 
erties. Due to storage constraints only a small fraction of the total number of 
pulses from PSR J0901—4046 were captured. 


Astrometry 


We compare the positions of radio sources in our combined images to those 
of the Rapid ASKAP Continuum Survey (RACS) [30]. We achieve this by 
running the PYBDSF [61] source finder on a joint MeerK AT image (made prior 
to model subtraction from all of the 2 second data), and the RACS image 
of the corresponding region. We forego the use of the RACS catalogue and 
run PYDSF directly on the RACS image in order to minimize any differences 
that may be introduced by the use of differing source finders. Furthermore, 
we select only sources that are compact in both catalogues (represented by 
a single point or Gaussian component) in order to minimize any resolution 
biases. The offsets between 254 MeerKAT sources and their matched RACS 
counterparts are shown in Supplementary Figure 1. The mean offsets in RA 
and Dec respectively are 0.90 (+ 0.07)” and —0.55 (+ 0.05)”, where the quoted 
uncertainties are the standard error in the mean. The standard deviations of 
the offsets in RA and Dec are 1.04” and 0.84” respectively. 

RACS’ positional accuracy is compared to that of the International Celes- 
tial Reference Frame (ICRF) version 3 [62]. Using RACS as an intermediate 
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step we find negligible (—0.052’’) systematic offset in RA, and a —0.95” sys- 
tematic offset in Dec (approximately 1/6th the size of the restoring beam) 
between the MeerKAT positions and the ICRF v3. 


Location in the P — P parameter space. 


Pulsars dissipate their rotational kinetic energy through electromagnetic radi- 
ation and a wind of relativistic particles leading to a gradual increase (called 
spin-down) of the rotation period as they age. A combination of the spin- 
period and spin-down or period derivative gives an estimate of a pulsar’s age 
and magnetic field strength and places it in the P — P parameter space which 
is used to compare relative populations of pulsars and to potentially trace 
their evolution (see Figure 1). Over the lifetime of the pulsar, it is thought 
to cross over a death line into a region called the pulsar graveyard [63] where 
it is no longer expected to emit in the radio. The rotational parameters of 
PSR J0901—4046 place it in the upper right region of the P — P parameter 
space, below two of the three death-lines in Figure 1. These death-lines are 
dependent on conditions needed for pair production in the pulsar magneto- 
sphere. The conventional emission model for pulsars assumes the existence of a 
vacuum gap above the polar cap of a neutron star. In order to sustain pair pro- 
duction, the potential difference across this gap must be sufficiently large and 
this is no longer possible beyond this death line. As a result, pair production 
and consequently, radio emission ceases. 

Equation 9 from CR93 [22], illustrated by the solid line in Figure 1, rep- 
resents a model for very curved or twisted magnetic field lines with curvature 
radius of the magnetic field line, r. ~ R = 10° cm comparable to the radius R 
of the neutron star. It also suggests a polar cap size much smaller than that 
of the pure dipole field case. PSR J0901—4046 lies well beyond this death-line 
implying that this model is not viable. 

The model for curvature radiation from the vacuum gap (represented by 
Equation 4 in [52] and illustrated by the dot-dashed line in Figure 1), uti- 
lizes a multipole magnetic field configuration and relativistic frame-dragging 
effects to achieve the necessary potential difference for continued pair produc- 
tion. PSR J0901—4046 is also located beyond this death-line implying that an 
alternative model is required to explain the radio emission. 

The space-charge-limited flow (SCLF) model death line for curvature radi- 
ation (calculated using Equation 9 in [52] and illustrated by the dashed line in 
Figure 1) assumes a multipole magnetic field and lies just below the location 
of PSR J0901—4046 indicating that this model is potentially feasible. Addi- 
tionally, it has been proposed [64] that different equations of state of a neutron 
star possibly affect the death-line, thereby leading to heavier radio pulsars 
surviving beyond the standard death line. 

For the simplest assumption of a dipolar magnetic field configuration, the 
minimum dipole magnetic field strength B at the surface of a canonical pulsar 
is, 
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which corresponds to 1.4 10!4 G for PSR J0901—4046. This value is similar to 
the bulk of the magnetar population in the P — P parameter space (calculated 
using this same standard formula). Assuming that the source was born with 
an initial period much shorter than its current period, a simplified estimation 
of the approximate, or characteristic age, T can be made as, 


P 
7 OP 8) 
which is equivalent to an age of 5.8 Myr for PSR J0901—4046. A magnetar 
with a spin period of ~ 10 s and a braking index n = 3, would evolve to the 
spin period of PSR J0901—4046 in ~ 5 Myr. We note that even if n were as 
high as 5, it still gives an age of ~ 2 Myr. 


Average Pulse Profiles and Pulse Modulation 


In Supplementary Figure 2 we show the average pulse profiles formed from the 
observations taken on 2021-04-02 when we observed with both the UHF and L- 
band separated by less than 2 hours to be able to compare the behaviour quasi- 
simultaneously between the two frequencies. The profiles are similar at the 
two frequencies with a smooth leading edge rising up to the peak of the profile 
and then a second component apparent on the trailing edge. The UHF profile 
is less smooth at the peak and has a more distinct trailing component which 
may be to do with the increased modulation (see below) in the UHF single 
pulses. Although the L-band data is from a single epoch, it is representative 
of the average pulse seen on other days. 

We measured the full-width at half maximum W5 9 of both the pulses by 
fitting five and seven von Mises functions (the von Mises distribution resem- 
bles a Gaussian distribution, but is cyclic) to the L- and UHF-band profiles 
respectively, using the PSRSALSA package [65]. This analytic, noise-free descrip- 
tion of the profile allows a robust width measurement of the profile shape. 
The uncertainties on the measurements are estimated by bootstrapping, i.e. 
by repeatedly adding Gaussian noise with the standard deviation measured 
from the off-pulse region of the profile, over multiple iterations. The quoted 
uncertainty is the standard deviation of the iterations. We find that the Wso 
widths of the two pulses are identical to within the errors and also to within 
0.3% of the pulse width indicating that between the two bands there is no evi- 
dence for radius-to-frequency mapping. It is interesting to note that the Wso 
point lies just about the point where the trailing feature is present indicating 
that too is apparently not evolving dramatically with frequency. 

As can be seen in Figure 1 there are a number of different types of pulse 
shapes that are seen from PSR J0901—4046. To try and quantify the variability 
we calculated the modulation index and standard deviation as a function of 
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pulse phase from a series of pulses (see [66] for the definition of these terms). 
We again used the UHF and L-band data from 2021-04-02 and the data set 
was the same as that used for the pulsar timing. We formed a fully-frequency 
averaged stack of single pulses for each of the 30 minute sessions, two each 
at the two frequencies, using the ephemeris given in Table 1. RFI mitigation 
was done before frequency averaging and was the same as described for the 
timing analysis above except that in some cases the pulses were corrupted too 
much by the RFI and these pulses were removed. We then combined the two 
30-minute stacks to form a single stack for each frequency, resulting in about 
45 pulses in each case. The single pulses were recorded with 65536 phase bins 
across a pulse period to give a time resolution of approximately 1.15 ms. The 
subsequent analysis was all undertaken using the PSRSALSA package. The off- 
pulse baseline was removed from each single pulse, and in some cases there 
was a slope in the baseline across the pulse and this was removed for doing 
a linear fit to the off-pulse region. The data were then gated to the section 
around the pulse profile as seen in Supplementary Figure 2. 

There are stark differences between the modulation properties at the two 
frequencies, which were observed less than 2 hours apart. We note that a similar 
analysis that used all the L-band pulses over the 6 epochs in 2021 gives a result 
very similar to the one obtained from this single day L-band data set. We 
only have this single epoch of UHF data, but given the stability of the L-band 
properties it is likely that this is representative of the UHF data and even if it 
is not, no similar outlier is seen in all the previous L-band observations, which 
still makes this difference very interesting. We note that the S/N at UHF is 
significantly higher than the L-band data and this is supported by the greater 
flux measured in the imaging data at UHF (weighted mean of 169.3 + 14 mJy 
beam~!, compared to 89.3 + 2.7 mJy beam~! at L-band). The modulation 
and standard deviation at L-band are reminiscent of many other pulsars, with 
the former peaking at the pulse edges and the latter somewhat tracing the 
intensity and neither being particularly high. This correlates well with the very 
low pulse jitter discussed in the timing section and confirms that despite the 
fact that the individual pulses are radically different, they form a stable average 
after a surprisingly small number of pulses. In contrast the UHF data show 
very strong modulation with parts of the leading edge exceeding a modulation 
index of one and a clear dip where the transition to the trailing component 
happens in the average pulse profile. We also see that the trailing component 
is much more highly modulated than it is at L-band. The standard deviation is 
also very high and traces the average profile indicating that the overall pulse- 
to-pulse flux is changing significantly. The imaging data also shows a greater 
standard deviation in flux between pulses at UHF (15 mJy) than at L-band 
(5 mJy) confirming the phase resolved analysis. 

Observations of the 23-second pulsar [19] at 350 MHz showed significant 
modulation in the pulse profile, that had not been seen at other frequencies, 
which was particularly evident in the leading component. However, these data 
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were not contemporaneous and the profile was different than that seen at 
higher and lower frequencies [67]. 


Polarimetric Calibration 


The complex channelized data for every pulse from PSR J0901—4046 were 
processed using a custom-made pipeline to generate beamformed, full-Stokes 
dynamic spectra (see “Voltage data capture pipeline” for more details). Then, 
we calibrated the polarization using the methodology given below. We assume 
that any leakage between the two polarization hands affects only Stokes V 
(defined as V=LL-RR using the PSR IEEE convention78). We also assume 
that the delay between the two polarization hands affects only Stokes Q and 
U. The calibration we apply ignores second-order effects. We perform a brute- 
force search for the rotation measure that maximizes the linear polarization 
fraction by using the PSRCHIVE tool rmfit [42]. With rmfit we select a range 
of rotation measures to search, in a number of trial steps. The delay between 
the polarization hands approximately manifests as an offset from the true 
rotation measure of the source, assuming the delay is frequency-independent. 
While this method provides a correct rotation measure, it still does not cali- 
brate the polarization position angle (PPA) to its absolute value. The complex 
gain solutions that are computed by the Science Data Processor pipeline, also 
account for the phase offset and delays between the two dipoles of the receiver 
and the phase offsets and delays with respect to the sky [68]. Since the data 
for each single pulse were only 1.5 seconds long, we ignored the phase offsets 
because of the rotation of the feed with respect to the sky. In order to confirm 
whether we have absolute PPA calibration, we also captured complex voltage 
data for a very well-known, bright pulsar PSR J0437—4715 whose polarization 
is very well measured at 1.4 GHz. We processed the data for PSR J0437-4715 
with same pipeline and corrected the full-Stokes data for rotation measure 
with the best known RM value taken from [69]. Then we compared the result- 
ing polarization of PSR J0437—4715 with a fully calibrated profile taken from 
the publicly available EPN pulsar database (http: //www.epta.eu.org/epndb/) 
[70]. Extended Data Figure 6 shows examples of fully calibrated single pulses 
of PSR J0901—4046. It is clear from the 1284 GHz pulses that the PPA follows 
a S-shaped sweep akin to PPAs observed for canonical radio pulsars. This is 
more obvious in the phased resolved PPA histogram shown in Supplementary 
Figure 3. 


Rotating Vector Model fits 


The Rotating Vector Model (RVM; [71]) provides a geometric interpretation of 
the variation of the PPA with rotational phase. To increase the sensitivity 11 
consecutive pulses were summed together. The resulting polarimetric profile 
is shown in Supplementary Figure 4. A rapid swing in PPA is shown near the 
peak of the profile, with the steepest part of the swing coinciding with the 
transition between the two main components of the profile. There is a distinct 
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dip in the degree of linear polarization. Although such a dip can in principle be 
the result of mixture of orthogonal polarization modes, the same orthogonal 
mode appears to dominate in this pulse longitude region for all the analyzed 
pulses (see the PPA distribution in Supplementary Figure 3). 

The PPA swing as observed in individual pulses is variable, especially where 
the degree of linear polarization is low. Nevertheless, the steep swing is a stable 
feature. The PPA swing of the 11 summed pulses was fitted with the RVM 
following the methodology described in [72]. The best fit (red line in the bottom 
panel of Supplementary Figure 4) is relatively poor (the reduced y? is 16), with 
significant deviations from what can be modelled with the RVM. The inability 
of the RVM to reproduce the observed PPA swing accurately means that 
the model parameters are poorly constrained. The magnetic inclination angle 
(with respect to the rotation axis) is unconstrained. What can be constrained 
is that the line of sight passes very close to the magnetic axis, with an impact 
parameter 8 S 0.2°. 

A small impact parameter is expected for a slowly rotating pulsar, as it will 
have an extremely large light cylinder, and a correspondingly small magnetic 
pole defined by open magnetic field lines on which particle acceleration and the 
production of radio emission is expected. The PPA swing therefore suggests 
that the observed radio emission originates close to the magnetic pole of the 
neutron star. 


Sub-pulse width evolution 


On top of the well-defined quasi-periodic behaviour and substructure shape, 
multiple pulses show a complex evolution in time. An example of such a pulse 
can be seen in the ’quasi-periodic’ panel of Figure 2. Here the pulse starts 
with relatively fast oscillations and the oscillations appear to slow down closer 
to the pulse peak. Interestingly, most quasi-periodic pulses show such ‘driven’ 
oscillation suggesting some sort of a ‘charge and discharge’ mechanism. Such 
behaviour may be too complex to be easily identifiable within a single ACF, 
which in this case is more sensitive to the shorter period. In order to investi- 
gate the evolution of the components and their quasi-periodicity we perform 
a continuous wavelet transform (CWT) on pulses detected with the APSUSE 
instrument as those data possess a finer time resolution (76 ys). This approach 
allows us to analyze our signal in both the time and frequency domains. The 
CWT of signal z(t) can be computed using the equation, 


+ inf 


CWT (a,b) = = Pte" (—) dt, (4) 


where ~(t) is the mother wavelet. The wavelet can then be stretched or com- 
pressed by a scale factor a and shifted by the translational parameter b. By 
varying the scale factor, the wavelet transform can be sensitive to features 
with different sizes. By applying the shift, the transform can localize the signal 
details in time. In our analysis we use two wavelets to extract different features 
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from the pulse profiles. The first order derivative of the Gaussian wavelet, later 
referred to as the G1 wavelet, at a time t given by 


We(t) = Cat exp (-t?) , (5) 


is used to detect sharp changes in the pulse profiles. By using a wavelet best 
suited to designing sharp transitions in the signal, we can obtain the start and 
end times of the quasi-periodic emission and derive the widths of individual 
components. To better understand the evolution of the quasi-period, we detect 
individual components and localize them in time. We can distinguish separate 
sub-pulses by using the Ricker wavelet of the form 


2 , -? 
bal) = 5 (1- Bex (>) (6) 

Supplementary Figure 5 shows an example wavelet transform analysis per- 
formed for one of the pulses with a strong quasi-periodic behaviour. The 
de-dispersed signal at a time resolution of 76 jus is first convolved with a win- 
dow of size 64 samples. This reduces the amount of unwanted noise and makes 
the structures of interest more prominent, without negatively impacting the 
time resolution. The convolved signal is then high-pass filtered to remove the 
overall pulse envelope. In this case, we use an ideal high-pass filter and fully 
filter-out frequencies below 5 Hz. The resulting convolved and filtered signal, 
shown in blue in all the panels of Supplementary Figure 5, is passed through 
the two wavelet transforms described above. In the case of the G1 wavelet the 
trough-to-peak transition results in a positive value of the wavelet transform 
at the point of the transition. Conversely, the peak-to-trough transition gives 
a negative value of the wavelet transform. When using the Ricker wavelet, the 
transform has its maximum at the centre of every peak. For both transforms 
we then use a simple peak finding algorithm to find and extract features of 
interest. The bottom row of Supplementary Figure 5 shows the resulting peak 
and trough width (left panel) and peak-to-peak period (right panel) measure- 
ments. For this pulse, the sub-pulse widths and separations are relatively stable 
and the averaged quasi-period 82.90 ms is consistent with the value of 79.7 ms 
obtained using the ACF method. 

We repeat this process for a sample of the detected pulses. Of those, 39 had 
at least 1 clearly defined sub-pulse that we were able to measure the width of. 
Supplementary Figure 6a shows the distribution of 145 sub-pulse widths that 
we have obtained with the wavelet method. The solid vertical line represents 
the median sub-pulse width of 49.00 ms. The dashed vertical lines represent 
one median absolute deviation limits. We can see that the majority of sub-pulse 
widths are concentrated around the median value, with only a small number 
of detections beyond 100 ms. The largest width detected within our sub-pulse 
sample was 120.43 ms, which is close to 50% of the FWHM at L-band reported 
in Table 1. 

Supplementary Figure 6b shows the relationship between the quasi-period 
of the sub-pulse oscillations and the width of the components of the individual 
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pulses. To reduce the complexity, only median values of component widths and 
their quasi-periods for every pulse are included. The solid black line represents 
a theoretical behaviour where the width of the sub-pulse is exactly half of the 
quasi-period. This describes a sub-pulse emission that spends half of the period 
in the ‘on’ state, followed by half of the period in the ‘off’ state. For our sample 
of pulses, only a handful can be classified as existing in such a state, with the 
divergent behaviour away from this half-period line for higher quasi-periods. 


Extant data 


We looked at archival synthesis images from various imaging surveys in the 
Southern Hemisphere. We searched for point source like emission in the closest 
fields in the TIFR-GMRT Sky Survey (TGSS) observed on 2016 March 15 at 
150 MHz [73], the Sydney University Molonglo Sky Survey (SUMSS) observed 
between February 2003 and July 2007 at 843 MHz [74] and the Rapid ASKAP 
Continuum Survey (RACS) observed on 2019 April 30 at 887.5 MHz [75]. 
We did not find any significant emission at the best known radio position 
of PSR J0901—4046 in all these surveys and report a 3-sigma upper limit of 
10.2 mJy, 6.6 mJy and 1.5 mJy for TGSS, SUMSS and RACS respectively. 
Taking the pulse-averaged flux density at 1.4 GHz as 200 Jy and the mea- 
sured spectral index of —1.7, i.e. approximately consistent with those of known 
rotation-powered pulsars [76], the expected flux densities would be 8.9 mJy, 
0.47 mJy and 0.43 mJy for TGSS, SUMSS and RACS, respectively. Addi- 
tionally, we examined the literature for pulsar time-domain surveys that had 
covered the PSR J0901—4046 field in the past. The field was indeed observed 
as part of the Parkes Multibeam Pulsar Survey (PMPS, [31]) in 1999 March 
(pointing centre offset by about 2.2 arcmin) and there were three further PMPS 
pointings near it in 1999 July/August, albeit offset by about 12 - 14 arcmin. We 
obtained the data and searched them for single pulses and periodic emission, 
both blindly and by folding using the known ephemeris. We repeated the anal- 
ysis without applying a high-pass filter, as was originally done in the default 
PMPS data processing, on the dedispersed time series. Neither the single-pulse, 
nor the periodicity search resulted in a detection. This is not surprising, as the 
data from the primary pointing are heavily affected by RFI. There was also 
a high-pass filter consisting of a 2-component RC time constant included in 
the data acquisition system which meant that longer timescale structure than 
0.9 seconds would be removed from the data. However, as the pulse width is 
less than this, some harmonics would still get through into the spectrum. The 
sensitivity penalty suffered because of this would therefore depend to a large 
extent on the pulse width, as well as the period. We estimate the 8-sigma flux 
density upper limit assuming a 1% duty cycle to be about 0.3 mJy. This too 
suggests that the non-detections in archival data (both synthesis imaging and 
time-domain) are not unexpected. 
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Energy distributions 


We examine the possibility of various types of pulses corresponding to varia- 
tions in the emission mechanism by computing pulse energy distributions. For 
this analysis we group the quasi-periodic and partially nulling pulses together, 
as all quasi-periodic pulses are partially nulling but not vice versa. It is not 
straightforward to separate these two classes without some sort of a qualitative 
metric which is beyond the scope of this paper. The average on-pulse energy 
distribution for each pulse archetype across both frequencies and all epochs 
is shown in Supplementary Figure 8. We use the PSRSALSA software suite to 
model pulse profile envelopes for the quasi-periodic and partially nulling pulses 
using a low-pass filter (see Supplementary Figure 7) to estimate the energy a 
pulse would have without the dropouts in power (Panel 2 from the top in Sup- 
plementary Figure 8). We lose ~ 40% of the energy to the observed dropouts. 
When this loss is accounted for by the profile envelopes, the energy distribu- 
tion is not very different from the distributions of the ‘normal’ or ‘split-peak’ 
pulses. These oscillations are likely linked to the radio emission production 
itself, or a periodic absorption mechanism that can suppress the radio emis- 
sion. The brightest observed pulses are classified under the split-peak category 
and are UHF-detections which is consistent with the observed spectral index 
reported in Table 1. 


Optical Follow-up 


High speed photometry of the 17th mag Gaia optical source near the radio 
position of PSR J0901—4046 was obtained on 7 and 9 January 2021, using 
the 1-m telescope of the South African Astronomical Observatory and the 
SHOC camera. On 7 January 2021, the source was observed with no filter and 
5 second continuous exposures for 2.66 hr (start time BJD 245 9222.41766). 
The second observation on 9 January 2021 again used no filter, and exposure 
times of 10 seconds continuous for 4.13 hr (start time BJD 245 9224.37231). 
The second observation showed some evidence for long term variations on time 
scales longer than the observation length, but no photometric modulations 
were seen at 76 seconds. 

Spectra of the same Gaia optical counterpart candidate was undertaken 
with SALT [77] on 2021 February 2. Two consecutive 1200 s exposures 
were obtained, beginning at 01:25:53 UTC, with the Robert Stobie Spec- 
trograph [78] which used the PG900 VPH grating, covering the region 
3920-6990 A at a mean resolution of 5.7 A with a 1”5 slit width. The spectra 
were reduced using the PYSALT package [79], (https: //astronomers.salt.ac.za/ 
software/pysalt-documentation/), which corrects bias, gain, amplifier cross- 
talk and cosmetic defects and finally mosaics the three CCDs comprising the 
detector. The spectral extraction, wavelength calibration and background sub- 
traction were all undertaken using standard IRAF routines, as was the relative 
flux calibration. The latter was achieved using an observation of EG21, taken 
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on 31 January 2021. The SALT spectrum of PSR J0901—4046 was flux cali- 
brated using an observation of the spectrophotometric flux standard, EG21, 
observed with the same grating setup two nights prior. Due to the inherent 
design of SALT, in particular the moving and variable entrance pupil, it is 
impossible to determine absolute fluxes of spectra. However, the relative fluxes 
and shape of spectra that are calibrated are reliable. In Supplementary Figure 
9 we show the spectrum of the Gaia optical counterpart candidate, which is 
typical of a reddened A-type star and therefore unlikely to be associated with 
PSR J0901—4046. 
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Supplementary Figure 1: Offsets between the MeerKAT positions and those 
of the Rapid ASKAP Continuum Survey [80]. We find negligible systematic 
offset of —0.052” in RA, and a —0.95” systematic offset in Dec between the 
MeerKAT positions and the ICRF v3. 
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Supplementary Figure 2: Pulse modulation properties of PSR J0901—4046 
at UHF and L-band. A comparison between the modulation (blue), standard 
deviation (orange) and average pulse profiles (dark) for the L-band (top panel) 
and UHF (bottom panel) data taken on 2021-04-02. There are some baseline 
variations still in the data due to the long period and so the standard deviation 
is truncated at a value of 0.15 and modulation indices are only plotted when 
the error is less than 0.1. The profiles are scaled to have a peak amplitude of 
1 and off-pulse mean of 0. 


Springer Nature 2021 TX template 


52 A radio-loud neutron star with an ultra-long spin period 


100 


-100 


Power (arb. units) 


Time (s) 
Supplementary Figure 3: Phase resolved polarization position angle his- 
togram. The top panel shows the phase resolved PPA histogram. The heat 
map of the PPA clearly traces out a S-shaped sweep reminiscent of the PPA 
swing seen in canonical radio pulsars. Bottom panel shows the average profile 
of the 23 pulses of PSR J0901—4046 added together. 
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Supplementary Figure 4: The polarimetric profile of PSR J0901—4046 
obtained by summing 11 consecutive individual pulses recorded on the first of 
February 2021. Top panel: The normalized Stokes J profile, the degree of linear 
polarization and Stokes V are shown as the black, red and green curves respec- 
tively. Bottom panel: The position angle (points with error bars) as detected 
with 20 confidence and the best fit RVM (red line). 
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Supplementary Figure 5: Example wavelet transform analysis for one of the 
pulses exhibiting a strong quasi-periodic behaviour. Two wavelets were used 
to estimate the emission properties. The first order derivative of the Gaussian 
wavelet was used to estimate the positions of peak-to-trough transitions, which 
provided an estimate of the individual sub-pulse widths and their evolution 
with time. The Ricker wavelet was used to determine the distance between 
consecutive sup-pulse peaks, providing independent and consistent estimates 
of quasi-periods to those obtained with the ACF method. 
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Supplementary Figure 6: Sub-pulse width distribution and the quasi-period 
— sub-pulse width relationship. Top panel: Sub-pulse width distributions for 
widths estimated using the wavelet method. The solid vertical line represents 
a median sub-pulse width of 49.00 ms. Bottom panel: Quasi-period — sub- 
pulse width relationship. The solid grey line represents a behaviour where 
the sub-pulses widths are half of their quasi-period, which is a case in e.g. 
sinusoidal oscillations. The dashed line represents a least squares fit, which 
shows a deviation from sinusoidal oscillations. 
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Supplementary Figure 7: Sample of fitted pulse profile envelopes for quasi- 
periodic and partially nulling pulses. The average energies of the profiles 
themselves as well as modelled envelopes are represented by the solid and 
dashed lines respectively. 
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Supplementary Figure 8: Pulse energy distribution of the different 
archetypes of pulses. ~ 40% of the energy is lost to the dropouts/dips seen 
in the quasi-periodic and partially nulling pulses in panel 1. When this is 
accounted for by modelling the pulse envelope (panel 2), the resulting distri- 
bution is similar to the other distributions in panels 3 and 4. This suggests 
that the pulses with dropouts/dips are not drastically brighter than the other 


types. 
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Supplementary Figure 9: SALT spectrum of coincident Gaia source, ruled 
out as the likely optical counterpart candidate. See text for details. 


